clear all
clear matrix
clear mata

set more off

local path = "...folder name"

cd "`path'"

local dlwcprogram = "acf_mata_steel.do"

do `dlwcprogram'

use steel_plant_pf.dta, clear

merge m:1 size using steel_k_price.dta /*Add capital price*/
drop _merge

gen cap_value = size*cp*10000

/*generate log variables*/
gen y = log(production)
gen l = log(labor)
gen k1 = log(size)
gen k2 = log(cap_value) 
gen logiron = log(tpigiron/1000)
gen logscrap = log(tscrap/1000) 
gen age = year - birth
gen m = logiron

gen dsecond = 0
replace dsecond = 1 if second > 0 & second < = 100

gen p = logscrap /*use scrap steel as proxy*/
gen k = k2 /*price adjusted capital*/

*trim outliers
drop if labor < 1 | labor > 500000

*---------------creating variables---------------------------------------------*
* higher order terms on inputs k l a p

forvalues ll = 0(1)5{
         local kmax = 5 - `ll'
		 forvalues kk = 0(1)`kmax'{
		      local amax = 5 - `ll' - `kk'
			  forvalues aa = 0(1)`amax'{
			       local mmax = 5 - `ll' - `kk' - `aa'
				   forvalues mm = 0(1)`mmax'{
				    local pmax = 5 - `ll' - `kk' - `mm' - `aa'
					forvalues pp = 0(1)`pmax'{
				    gen poly`ll'`kk'`mm'`aa'`pp' = (l^`ll')*(k^`kk')*(m^`mm')*(age^`aa')*(p^`pp')
				     }
				   }
			  }
		 }
}


*------FIRST STAGE USING EXP AS INPUT  ----------------------------------------*
local regphi = "poly* dsecond second i.t i.owner i.region"

xi: reg y `regphi'
predict phi
predict epsilon, res
label var phi "phi_it 
label var epsilon "measurement error first stage

local regx = "l k m age i.t i.owner i.region"

xi: reg y `regx' if !missing(p)
for any l k m: gen OLSX=_b[X]
gen OLSConst=_b[_c]	
gen OLSa = _b[age]

*------------------------------------------------------------------------------*
xtset id time, monthly
gen phi_lag=L.phi
gen l_lag=L.l
gen k_lag=L.k
gen age_lag = L.age
gen m_lag = L.m
gen dsecond_lag = L.dsecond
gen second_lag = L.second
sort id t
gen const=1
drop if y==.
drop if l_lag==.
drop if l == .
drop if k==.
drop if m ==.
drop if m_lag == .
drop if phi==.
drop if phi_lag==.
drop if age ==.

*try ols initial value first

for any l k m: gen initialX=OLSX
gen initiala = OLSa
gen initialx = sqrt(initiall/(1-initiall))
gen initialz = sqrt(initialk/(1-initialk))

dlwc

matrix result = (beta_dlw[1,1]^2/(1+beta_dlw[1,1]^2),beta_dlw[1,2]^2/(1+beta_dlw[1,2]^2),beta_dlw[1,3],beta_dlw[1,4],beta_dlw[1,5],initiall[1],initialk[1],initialm[1])

putexcel A1 = mat(result) using acf_steel_lkm_kprice_link, sheet("iron_scrap") modify

local n = 2
local iter = 101

set seed 12345 

while `n' < = `iter'{

gen u = (1.2 - 0.8)*runiform() + 0.8 /*search around 80% to 120% of ols coefficients*/

forvalues i = 1(1)4{
gen u`i' = u[`i']
}

replace initiall = OLSl*u1
replace initialx = sqrt(initiall/(1-initiall))
replace initialk = OLSk*u2
replace initialz = sqrt(initialk/(1-initialk))
replace initialm = OLSm*u3
replace initiala = OLSa*u4

dlwc

matrix result = (beta_dlw[1,1]^2/(1+beta_dlw[1,1]^2),beta_dlw[1,2]^2/(1+beta_dlw[1,2]^2),beta_dlw[1,3],beta_dlw[1,4],beta_dlw[1,5],initiall[1],initialk[1],initialm[1])


putexcel A`n' = mat(result) using acf_steel_kprice_link, sheet("iron_scrap") modify

matrix drop result
drop u u1-u4

local n = `n' + 1
}
